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Abstract 

Spectral methods for compressible flows are introduced in relation to 
finite difference and finite element techniques within the framework of the 
method of weighted residuals* Current spectral collocation methods are put in 
historical context. The basic concepts of Fourier spectral collocation 
methods are provided. Filtering strategies for shock-capturing approaches are 
also presented. Fourier shock capturing techniques are evaluated using a one- 
dimensional, periodic as t rophysical "nozzle" problem. 
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1% Introduction 


Spectral methods may be viewed as an extreme development of the class of 
discretization schemes known by the generic name of the method of weighted 
residuals (MWR).^ The key elements of the MWR are the trial functions (also 
called the expansion or approximating functions) and the test functions (also 
known as weight functions)# The trial functions are used as the basis 
functions for a truncated series expansion of the solution, which, when 
substituted into the differential equation, produces the residual. The test 
functions are used to enforce the minimization of the residual# 

The choice of trial functions is what distinguishes the spectral methods 
from the finite element and finite difference methods. The trial functions 
for spectral methods are infinitely differentiable global functions. 
(Typically they are tensor products of the eigenfunctions of singular 
S turm-Liouvi 1 le problems.) In the case of finite element methods, the domain 
is divided into small elements, and a trial function is specified in each 
element. The trial functions are thus local in character, and well-suited for 
handling complex geometries. The finite difference trial functions are 
likewise local. 

The choice of test function distinguishes between Galerkin and 

collocation approaches. In the Galerkin approach, the test functions are the 
same as the trial functions, whereas in the collocation approach the test 
functions are translated Dirac delta functions. In other words, the Galerkin 
approach is equivalent to a least square approximation, whereas the 
collocation approach requires the differential equations to be satisfied 
exactly at the collocation points. 

The collocation approach is the simplest of the MWR, and appears to have 
been first used by Slater^ in his study of electronic energy bands in 
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metals, A few years later, Barta^ applied this method to the problem of 
torsion of a square prism. Frazer, et al.^ developed it as a general method 
for solving ordinary differential equations. They used a variety of trial 
functions and an arbitrary distribution of collocation points. The work of 
Lanczos^ established for the first time that a proper choice of the trial 
functions and the distribution of collocation points is crucial to the 
accuracy of the solution. Perhaps he should be credited with laying down the 
foundation of the orthogonal collocation method. This method was revived by 
Clenshaw, Clenshaw and Norton, and Wright. These studies involve 

application of Chebyshev polynomial expansions to initial value problems. 

Villadsen and Stewart^ developed this method for boundary value problems. 

The earliest applications of the spectral collocation method to initial 
value problems in partial differential equations were those of Kreiss and 
Oliger,^ (who called it the Fourier method) and Orszag,^ (who termed it 

pseudospec t ral) • Details can be found in Ref. 12. 

Spectral methods have been used on one-dimensional, compressible flow 
problems with piecewise linear solutions by GottlLeb, et al. J and by Taylor, 
et al.^. These reports demonstrated that spectral methods, when combined 

with appropriate filtering techniques, can capture one-dimensional shock waves 
in otherwise featureless flows. A different sort of demonstration was 

provided by Zang and Hussaini.^ They exhibited spectral solutions of 
compressible flows with non-trivial structure in the smooth regions. 

A systematic assessment of the accuracy of spectral methods for 
compressible flows is needed, however, to determine whether they are useful 
rather than merely feasible. Some comparisons with finite difference 

solutions were presented in Ref. 15, but the length constraint on that brief 
report precluded any detailed comparison. Cornille^ presented both spectral 


3 


and finite difference results for a step function solution of the inviscid 
Burgers' equation. His comparisons are not entirely satisfactory because the 
time discretization errors were not assessed and the finite difference results 
would surely have been better had the calculation employed a uniform grid. 

The purpose of the present paper is to assess the accuracy of Fourier 
spectral methods when used to capture shock waves. In particular, emphasis 
will be placed on the role of filtering which must be used for stability 
purposes and for the elimination of oscillations. One of the test problems 
used by Zang and Hussaini^ provides a non-trivial test of the filtering 
strategies. 

The restriction here to Fourier methods limits the discussion to periodic 
problems. Such methods are appropriate, however, for an important class of 
problems represented by flow around an aerofoil. In a companion paper^ , we 
address two-dimensional non-periodic problems and the use of Chebyshev methods 
in conjunction with shock fitting. 


II. Spectral Methods for Shock-Capturing 

Thus far the spectral collocation method has been the only type of 
spectral method applied to compressible flow problems. The present discussion 
will be confined to spectral collocation methods, with all future references 
to spectral methods implicitly referring to this specific type. Since 
spectral methods are a novel approach to aerodynamic flow computations, a 
basic introduction to their properties and implementation will be presented 
first. 
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Basic Fourier Spectral Concepts 

The potential accuracy of spectral methods derives from their use of 
suitable high-order interpolation formulae for approximating derivatives. 
Their efficiency has often depended on the use of Fast Fourier Transform 
techniques. An elementary example is provided by the model problem 

u t + u x = °* (1) 

with periodic boundary conditions on [0,2 tt] and the initial condition 

u(x,0) = sin(TT cos x) . (2) 


The standard collocation points are 


x . 
J 



J = 0 ,!,••• ,N-1 • 


( 3 ) 


Let uj denote the approximation to u(xj), where the time dependence has 
been suppressed. Then the spatial discretization of Eq. (1) is 



/v 

9 u 

3 x j 


(4) 


where the right-hand-side is determined as follows. First, compute the 
discrete Fourier coefficients 


u k 


1 N" 1 -ikx 

1 “je 
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(5) 


Then the interpolating function 
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N/2-1 

u(x) = l G k e ikx (6) 

k= -N/2 

can be differentiated analytically to obtain 

N/2 - 1 - ifcx, 

l ik u k e J . (7) 

k= -N/2+1 

(The term involving k = -N/2 makes a purely imaginary contribution to the 
sum and hence has been dropped.) Note that each derivative approximation uses 
all available information about the function values. The sums in Eqs. (5) and 
(7) can be obtained in 0(N £n N) operations via the Fast Fourier Transform 
(FFT) . 

An illustration of the superior accuracy available from the spectral 

method for this problem is provided in Table I. Shown there are the maximum 
errors at t = 1 for the spectral method as well as for second-order and 

fourth-order finite difference methods. The time discretization was the 
classical fourth-order Runge-Kutta method. In all cases the time-step was 
chosen so small that the temporal discretization error was negligible. 
Because the solution is infinitely smooth, the convergence of the spectral 

method on this problem is more rapid than any finite power of 1/N. (The 
error for the N = 64 spectral result is so small that it is swamped by the 
round-off error of these single precision CDC Cyber 175 calculations.) This 
type of convergence is usually referred to as exponential. In practical 

applications the benefit of the spectral method is not the extraordinary 
accuracy available for large N but rather the small size of N necessary 

for a moderately accurate solution. 


3u 
<Tx j 
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Of course, full-scale fluid dynamical problems are considerably more 
complicated than those described by Eq. (1). Nevertheless, for quite a few 

IQ IQ # 

three-dimensional, viscous, incompressible problems these complications 

have not prevented spectral methods from displaying the sort of accuracy and 
economy suggested by the results in Table I. 


Table !• Maximum Error for a 1-D Periodic Problem 


N 

Fourier 

Spectral 

2nd-0rder 

Finite 

Difference 

4 th-0rder 
Finite 
Difference 

8 

1.62 (-1) 

1.11 (0) 

9.62 (-1) 

16 

4.97 (-4) 

6.13 (-1) 

2.36 (-1) 

32 

1.03 (-11) 

1.99 (-1) 

2.67 (-2) 

64 

9.55 (-12) 

5.42 (-2) 

1.85 (-3) 

128 


1.37 (-2) 

1.18 (-4) 


The complications to Eq. (1) relevant to the Euler equations are: (1) 

non-constant coefficients, (2) non-periodic boundary conditions, (3) a severe 
explicit time-step restriction, and (4) nonlinearities causing shocks. The 
first three of these are also present for the incompressible Navier-Stokes 
equations. Non-constant coefficient problems may be technically unstable but 
this can be controlled by filtering techniques. ^ Highly accurate solutions 
to problems with non-periodic boundary conditions can be obtained by using 
Chebyshev polynomials in place of trigonometric functions. (The use of 
Chebyshev polynomials is discussed in the companion paper.^) A widely 
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applicable technique for surmounting the explicit time-step restriction has 
yet to appear. However, since this difficulty is algorithmic rather than 
conceptual, it should be resolved eventually. Finally, the most serious 
hurdle is surely the global oscillations arising from the presence of shocks 
in the interior of the computational domain. This also seems to be a problem 
in filtering so an extended discussion of filtering in this context is 
provided in the next subsection. 

Filtering for Fourier Spectral Methods 

Several types of filtering operations are employed in spectral methods: 

1) Pre-processing . The initial condition is filtered, usually in Fourier 
space. New grid point values are obtained from 

N/2-1 * ikx . 

Uj = l a(277k/N) u k e J , (8) 

k=-N/2 

where a(0) is a non-negative function defined on [ — tt , tt ] such that 

a(0) - 1 for 1 0 1 — 0 and o(0) ** 0 as 1 0 | + tt. The coefficients u^ 

could be either the discrete Fourier coefficients of the original initial 

condition u(x,0) or its continuous Fourier coefficients defined by the usual 

integral. 

2) Derivative filtering . In the computation of spatial derivatives the 
term ik in Eq. (7) is replaced with ika (2TTk/N) . 

3) Solution smoothing . At regular intervals in the course of advancing 
the solution in time, the current solution values are smoothed in Fourier 
space in a manner described by Eq. (8). 

4) Cosmetic post-processing . If only weak filtering is needed in order 
to stablize the computations, it may be necessary to smooth the solution as 
above for display purposes only. 
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5) Artificial viscosity . Just as with finite difference computations, 
the equations can be modified with either a linear or non-linear artificial 
viscosity. 

The theoretical reasons for employing pre-processing are confined to 
linear problems with discontinuous solutions. 2 ^ None of the calculations 
reported here employed pre-processing. 

Derivative filtering has been proposed 20 as an inexpensive way to deal 
with the instability of Fourier spectral approximations to non-constant 
coefficient problems. For a (linear) constant coefficient equation, 
derivative filtering is identical to a similar pre-processing and to solution 
smoothing. (An alternative but more costly approach is to re-write the 

r\ r\ 

equations in skew-symmetric form. z ) Our own experience with derivative 
filtering on more complicated problems has been unfavorable. Its absence has 
not led to any noticeable instability and its presence has often caused 
troublesome oscillations in the solution. All the results in this paper have 
been obtained without any sort of derivative filtering. 

To date, the most successful smoothing strategy to control the 
oscillations which develop when shock capturing methods are used has been 
solution smoothing. With this approach, however, one must choose not only the 
particular filter to be used, but also the frequency of application. As with 
artificial viscosity methods for finite difference techniques, fine tuning of 
the smoothing process must be done in order to obtain the best results. 

Cosmetic post-processing was introduced by Gottlieb et al.^ who used it 
on problems where only weak smoothing was needed for stability. A strong 
filter is then used only to make the solution presentable, not to modify the 
solution. 
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Finally, an artificial viscosity can be applied not only to stabilize the 
computations but to smooth them out as well# These can be viewed as non- 
linear filters. They are not, however, well matched to spectral computations 
since they introduce a finite order truncation error which negates the use of 
the high order interpolation used by the spectral method. 

Some of the smoothing functions that have been employed in our 
calculations are 

o* ( 0 ) = — (Lanczos) (9) 
a(6) = (1 + cos 0) (raised cosine) (10) 
o(0) = aj(35 - 84 a Q + 70a^ - 20aJ) (11) 

(sharpened raised cosine) 


where 


is the raised cosine given by Eq 


o(0) 


1 

( 0 “ 7T ^4 

c 


( 10 ), 

M < 9 C 

e c < |0| < tt 


(quartic taper) 


o(0) 


1 l«l < » c 

-«<M - e.) 4 


( 12 ) 


( 13 ) 


e 


0 c < |0 | < TT 

(exponential cut-off) 
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These filters are listed above in order of decreasing strength and are 
illustrated in Fig. 1. The Lanczos and raised cosine filters are classical. 
Equation (11) represents one of a number of standard formulae for sharpening a 
basic filter. The exponential cut-off has been proposed specifically as a 
filter for use in spectral methods. 21*13 ^ popular choice for the cut-off 

function has been the exponential cut-off. Usually 0^ lies between tt/2 

and 5 tt / 6 and a is chosen so that o(tt) is 0(10“^) or smaller. Some 

2 2 

successful results have been reported for this filter on linear problems. 

The choice of filter will determine which Fourier frequencies will be 
modified. The k = 0 component of the Fourier decomposition is the only one 
which contributes to the average value. Thus, in order to preserve the mean 
value of the solution to a conservation law, o(0) =1 is required. The 
filters shown in Fig. 1 all have this property. The effect of the filter on 
the high frequencies, however, is usually more difficult to assess in 
nonlinear problems • 

The raised cosine (also known as the von Hann window) admits a simple 
physical interpretation. It is algebraically equivalent to 


u . , + 2u . + u . , . 

u. ♦ -a i — J — ill 

J 4 


. „ . (ax ) 2 “j-i - 2 "j * “j+i 

j 4 (Ax) 2 


(14) 


This is clearly a second-order artificial viscosity term with the coefficient 


(Ax) 2 

4(At/N f ) ’ 


( 15 ) 


11 


where Nf is the number of time-steps between applications of the filter. 
Figure 1 suggests that the other types of filters amount to non-physical 
viscosities: they damp preferentially the high frequency components of the 
solution but in a different manner than a physical viscosity* 

The frequency of applying the filter is analogous to the selection of the 
size of the artificial viscosity for finite difference methods. Applied too 
often a strong filter like the Lanczos filter will unacceptably smear out a 
shock. Frequent applications of a weak filter such as the exponential cut-off 
on the other hand may not be enough to stabilize a solution let alone 
eliminate oscillations. At present it is still necessary to determine the 
best filter and its frequency of application by trial and error. 

Ideally a filter should suppress the oscillations arising from a 
discontinuity while retaining the accuracy of the spectral method in the 
structured regions of the flow. None of the filters described above is 
entirely satisfactory. A simple illustration is furnished in Figure 2. Shown 
there are several approximations to one period of a function composed of a 
square wave plus a small, narrow Gaussian. The approximations are based on 
the interpolating polynomials generated by Eqs. (5) and (6) with N = 32. For 
the bottom two parts of Figure 2 the Fourier coefficients have been modified 
by the raised cosine and the exponential (with 9 c = 2 tt/ 3) filters, 
respectively. The oscillations of the interpolating polynomials between the 
grid points are significant for real problems because they are indicative of 
the dangerous high-frequency modes that will interact with variable- 
coefficient and non-linear terms. 

The familiar dilemma illustrated by Figure 2 led to the adoption of non- 
linear filters in central difference schemes. As with artificial viscosities, 
the non-linear filters currently available introduce a truncation error of 
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finite order into the overall scheme, apparently precluding them from 
converging exponentially fast even in the smooth regions of the flow. These 
considerations imply that spectral shock-capturing solutions of compressible 
flows must be examined critically. 


III. Results for Fourier Spectral Shock-Capturing 

Discussions of shock-capturing techniques are easiest for the Fourier 
spectral methods. The discrete operator is simple, the collocation points are 
uniformly distributed and the boundary conditions pose no difficulty. 
However, there does not appear to be any interesting one-dimensional 
aerodynamic problem with periodic boundary conditions. One non-trivial test 
case for spectral shock-capturing techniques uses an approximate set of 
equations derived by Woodward 2 ^ for studying the time development of shock 
waves in a spiral galaxy. The equations which describe an isothermal gaseous 
component in a very thin disk galaxy are 

u + (pu) = 0 

t <f> 

(pu) t + [p(u 2 + a 2 )]^ = 20(v-v Q )p + pe sin <f> (16) 

2 

(pv) t + (Puv)^ = - (^)(u - u Q )p. 

The boundary conditions are periodic in (j> with period 2 tt. The last term 
in the middle equation represents the gravitational forcing of the gas from 
the spiral field of the much more massive stellar component. These equations 
have some significant differences from the Euler equations, notably the 
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forcing term and the asymmetrical role of the velocity components. A more 
detailed explanation of this physical problem and the approximations used in 
deriving Eq. (16) is available in Ref. 24. 

The parameter e in Eq. (16) is a dimensionless measure of the strength 
of the gravitational forcing. If this forcing is sufficiently strong, then 
the steady-state solution to this astrophysical problem contains a shock. 
Behind the shock there is a region of rapid decompression and further 
downstream occurs some features due to the second harmonic of the forcing 
term. The steady-state solution, then, is more complex than some standard 
test problems (such as the one-dimensional shock tube) whose solutions are 
merely piecewise linear functions. The challenge for the spectral method is 
to capture the shock and to suppress its attendant oscillations without also 
destroying the remaining structure of the solution. 

The specific test problem uses (in units which are not of interest here) 
a = 8.56, ft = 21.37, v Q = 115, e = 72.92, k = 26.75, and u Q = 13.5. 

The spectral calculations use the steady-state solution to Eq. (16) as 
the initial condition. The transients that would be generated by some other 
initial condition take a very long time to damp out because spectral methods 
have very low inherent damping and there are no boundaries out of which 
transients can convect. The temporal discretization uses a second-order 
Adams-Bashforth predictor followed by a third-order Adams-Moulton corrector. 
It has a CFL limit of 3 / ( 2 tt ) • The calculations used a CFL number of 0.2 and 
were run for 500 time-steps. We seek to illustrate the performance of 
filtering procedures ranging from the very weak to the very strong. The 
specific examples to follow are representative of many dozens of calculations 
in which the type of filtering, the form of the filtering function, and the 
frequency of filtering were each varied over a wide range. 
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Figure 3 shows the effect of applying the weak exponential cut-off filter 
(with 0 C = 0.7 tt and a = 5) every 50 time-steps. Both high- and low- 

frequency oscillations are quite evident. If a Lanczos cosmetic 
post-processing step is added, then the results shown in Fig. 4 are 
obtained. As one would expect, the high-frequency oscillations are under 
control but the low-frequency errors remain substantial. The failure of the 
weak filtering strategy on this problem (even when applied every time step) 
may appear puzzling in view of the success a similar strategy has achieved on 
a standard shock tube problem.^ We attribute the difference to the strong 

expansion in the post-shock region in the present problem. Without 
appreciably stronger smoothing something resembling an expansion shock will 
form in this region. The dip in the density plot is not smoothed by a weak 
filter and it grows until the density becomes negative and the calculation 
breaks down. On the other hand, if the forcing parameter £ is reduced to 
31.0 then the weak exponential filter is sufficient to produce a stable 
computation. However, the shock here is quite weak (the density ratio is 
1.38). The more difficult problem for £ = 72.92 evidently requires more 
drastic filtering. 

Stable computations are obtained by applying the raised cosine filter 
every 50 time-steps. As Fig. 5 indicates, the expansion is now adequately 
controlled and the shock is captured in two points. But notice that there are 
some low-frequency errors in the vicinity of the second harmonic near a spiral 
phase of 270°. 

An alternative form of strong filtering is to add a nonlinear artificial 
viscosity. In his second-order MacCormack's method calculations for the two- 
dimensional version of this problem Liebovitch^ used a fourth-order 
viscosity. In the density equation this was proportional to 
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where 


v 


i+ 1/2 




v 


i- V 2 



P i-1 


), 


V 


i+ V 2 




(17) 


(18) 


Similar terms were used for the two momentum equations with the appropriate 
momentum variable replacing the density variable in Eq. (17)* Figure 5 shows 
the results of employing this nonlinear artificial viscosity instead of any 
type of linear filtering. Compared with Fig. 5 there is now one more point in 
the shock but the solution is smoother, particularly near the second 
harmonic. A cosmetic post-processing step which is one-sided near the shock 
would certainly sharpen the latter solution. 

The need for strong filtering to keep the computations stable limits the 
accuracy that is obtained for this problem. Table II shows the discrete L 2 
errors in the density excluding the shock region for both the strongly 
filtered and artificial viscosity solutions on grids of 16, 32 and 64 
points. Away from the shock both approaches give very similar errors and both 
show first order error decay rates. Nothing approaches "spectral accuracy", 
i.e., exponential convergence, has been obtained for this problem. 


Table 11% Comparison of L 2 Errors for Filtering and Artificial Viscosity 


N 

Filtering 

Artificial Viscosity 

16 

0.087 

0.098 

32 

0.036 

0.037 

64 

0.016 

0.017 
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A comparison with the preliminary results^ for this problem is in 
order. The earlier spectral result is given here in Fig. 7. The spectral 
method used a second-order linear artificial viscosity applied every time-step 
with the raised cosine solution smoothing added every 100 time-steps in 
physical space (see Eq. (14)) but modified for one-sided linear extrapolation 
for the two points straddling the shock. This latter feature accounts for the 
apparent sharpness of the shock. Unfortunately, that type of one-sided 
averaging was not conservative. Inspection of Fig. 7 reveals that the total 
mass has increased. Conservative versions of this one-sided averaging have 
not been as successful in producing sharp shocks. 

All these spectral results may be compared with the second-order 
MacCormack's solution displayed in Fig. 8. A nonlinear artificial viscosity 
as given by Eq . (17) was included in these calculations. Grid refinement 
studies for this method also reveal essentially first-order error decay. Away 
from the shock, this straightforward finite difference result appears to be 
superior to all of the spectral results. Upwind difference schemes in 
resolving the shock can do even better than MacCormack's method. Results for 
a second-order flux splitting method applied to this problem are given in 
Refs. 26 and 27. In terms of the error norm used in Ref. 26, the error in the 
density is 1.6% and 2.4% for the present spectral filtering and artificial 
viscosity calculations. For Van Leer's second order scheme it is only 0.6%. 


IV. Conclusions 

The astrophysical problem is the most challenging one-dimensional 
compressible flow problem for which spectral shock-capturing results have been 
reported. Cornille^ examined a scalar problem for which the MacCormack 
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finite difference method had sufficient implicit dissipation* In contrast, 
explicit dissipation had to be added to this finite difference method for it 
to handle the astrophysical problem* Cornille's spectral methods evidently 
required no explicit filtering on his problem because the shock is relatively 
weak and because they have an implicit dissipation similar to that of 

MacCormack's method. The shock tube problem of Gottlieb, et al.^ required 
only a weak exponential filter, whereas much stronger linear filtering was 
required for the astrophysical problem. The numerical examples of Taylor, et 
al.^ were simpler than the shock tube problem, consisting of only a single 

wave. The piecewise linear nature of these other test problems precludes 
their use for assessing the accuracy of a numerical method. 

The performance of the spectral method on the astrophysical problem 
should serve as a caution to those considering this approach. The only 
presentable solutions we obtained employed such drastic filtering that the 
accuracy of the method deteriorated to first-order. In a direct comparison 
finite difference methods produced better (and cheaper) solutions. 

We are not advocating that spectral methods be abandoned for compressible 
flows. There are circumstances under which they are quite competitive 

O Q 

indeed. Streett^ 0 has developed a spectral shock-capturing method for the 

two-dimensional full potential equation. His spectral method requires far 
fewer grid points than the best finite difference methods to achieve 
engineering accuracy for supercritical, lifting airfoils. Moreover, the 
application of spectral multigrid techniques^ makes the cost of this method 
comparable to that of finite difference methods. Streett's results are so 

accurate, albeit not exponentially so, because only a weak filter is 
required. The Gibbs oscillations are less severe for this problem than for 
the astrophysical problem because the solution itself is continuous at the 
shock; it is the derivative which is discontinuous. 
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Exponentially convergent spectral solutions to the Euler equations can be 
obtained by resorting to shock-f itting techniques. These methods are 
described for two-dimensional flows in the sequel to the present paper.^ 
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Figure Captions 


. Filters used in the spectral calculations. a. Raised Cosine; b. 
Lanczos; c. Sharpened raised cosine; d. Exponential cut-off. 

• Effect of several linear filters upon a periodic square wave plus a 
Gaussian bump. The dashed line indicates the exact function, the 
circles indicate the value of the pos t-f iltered approximation of the 
collocation points and the solid line indicates the interpolating 
polynomial. 

. Spectral solution to the astrophysical model (circles) computed with 
the exponential cut-off filter applied every 50 steps. The solid 
line is the exact solution. 

. Solution as in Fig. 2 but with a Lanczos cosmetic filter also 

applied. 

. Spectral solution to the astrophysical model with a raised cosine 
filter applied every 50 steps. 

» Spectral solution to the astrophysical model with a nonlinear 
artificial viscosity. 

* Spectral-finite difference solution using a weak artificial 

viscosity and raised cosine filter. One-sided linear extrapolation 
is used at the shock when smoothing is applied. 

► MacCormack solution to the astrophysical model. 
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